% This function is part of the replication package of "The Cross Section of Foreign Currency Risk Premia 
% and Consumption Growth Risk: Comment" by Craig Burnside, published on American
% Economic Review, vol.101, n.7, December 2011. The entire replication
% package is available on https://www.aeaweb.org/articles?id=10.1257/aer.101.7.3456

% 

% COPYRIGHT 2011 American Economic Association

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
% ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
% DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
% ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
% (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
% LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
% ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.


function S=hac_var(u,K1,K2)

% compute a HAC or VARHAC estimator of the long-run covariance matrix of GMM errors 

% INPUTS:
% u:  error matrix, T x n
% K1: maximum lag order considered for own lags
% K2: maximum lag order for other lags
% if K1=K2=0, HAC standard errors are computed

% OUTPUT:
% S: the estimated spectral density matrix of the errors

%#ok<*MINV>
%#ok<*AGROW>
%#ok<*WNOFF>
%#ok<*WNON>

if K1<K2 
    disp('Error: K1 cannot be less than K2')
    S=[] ;
else
    K=max([ K1 ; K2 ]) ;
    [T,n]=size(u) ;
    su=std(u,1) ;
    u=u./(ones(T,1)*su) ; % normalize the std. dev. of the GMM errors to 1 (undo this later)
    if K~=0
        v=u(K+1:T,:) ;
        for j=1:n 
            % work on the jth equation
            y=u(K+1:T,j) ;
            BIC=log(y'*y) ; % initialize the BIC
            A(j,:)=zeros(1,K*n) ; % initialize the jth row of the A matrix 
            for k=0:K1 
                % this is the number of lags of variable j to enter
                if k==0
                    x1=[] ;
                    minit=1 ;
                elseif k==1
                    x1=u(K:T-1,j) ;
                    minit=0 ;
                else
                    x1=[ x1 u(K+1-k:T-k,j) ] ;
                    minit=0 ;
                end
                for m=minit:K
                    % this is the number of lags of the other variables to enter
                    if m==0
                        x2=[] ;
                    elseif m==1
                        x2=u(K:T-1,[(1:j-1) (j+1:n)]) ;
                    else
                        x2=[ x2 u(K+1-m:T-m,[(1:j-1) (j+1:n)]) ] ;
                    end
                    x=[ x1 x2 ] ;
                    if m<=K2 
                        warning off ; 
                        b=(x'*x)\(x'*y) ;
                        warning on ; 
                        vv=y-x*b ;
                        bic=log(vv'*vv)+(k+m*(n-1))*log(T-K)/(T-K) ;
                        if bic<BIC
                            BIC=bic ;
                            sizb=max([k ; m])*n ;
                            bb=zeros(sizb,1) ;
                            if k>0
                                bb(j:n:(k-1)*n+j,1)=b(1:k)' ;
                            end
                            if m>0
                                for p=1:m
                                    bb((p-1)*n+1:(p-1)*n+j-1,1)=b(k+(n-1)*(p-1)+1:k+(n-1)*(p-1)+j-1,1) ;
                                    bb((p-1)*n+j+1:p*n,1)=b(k+(n-1)*(p-1)+j:k+(n-1)*p,1) ;
                                end
                            end
                            A(j,1:sizb)=bb' ;              
                            v(:,j)=vv ;
                        end
                    end
                end
            end
        end
        B=[ eye(n) (-A) ] ;
        Sv=(v'*v)/(T-K) ; 
        sumB=B(:,1:n) ;
        for j=2:(K+1) 
            sumB=sumB+B(:,(j-1)*n+1:j*n) ;
        end
        S=inv(sumB)*Sv*inv(sumB)' ; 
    else
        S=(u'*u)/T ; % no correction for serial correlation
    end
    S=S.*((su')*su) ; % scale the covariance matrix back up
end